function verify2(num_fig)
close all

if nargin == 0
    num_fig = 0;
end


global A delta

% define the initial condition and delta
u0 = [1; 0];
delta = 1e-2;
T = 1;

% define the matrix
A = @(tau) [    
    -tau.^3 ./ (1 + tau.^3), tau.^2 ./ (1 + tau.^3);
    tau ./ (1 + tau.^3), -1 ./ (1 + tau.^3)
    ];

% A = @(tau) [
%     exp(-tau) - 1/2 1;
%     0 -1];

c_1i = adj_evec(1, 0)' * u0;
c_2i = adj_evec(2, 0)' * u0;

[evec(1,0) evec(2,0)]
[adj_evec(1,0) adj_evec(2,0)]

% numerical solution
opts = odeset('abstol',1e-8, 'reltol',1e-8);
num = ode15s(@(t, X) A(delta * t) * X, [0, T / delta], u0, opts);
fprintf('num soln computed\n');

if num_fig
    subplot(2,1,1); plot(num.x, num.y(1,:));
    subplot(2,1,2); plot(num.x, num.y(2,:));
    pause;
    close all
end

% most basic solution
c1 = @(tau) c_1i * exp(beta(1, tau));
c2 = @(tau) c_2i * exp(beta(2, tau));

% correction

corr1 = @(tau) -c_2i * gamma(1, 2, tau) ./ (lambda(1, tau) - lambda(2, tau)) .* exp(beta(2, tau)) + ...
    ( c_1i * quad(@(rho) gamma(1,2, rho) .* gamma(2,1,rho) ./ (lambda(1, rho) - lambda(2,rho)), 0, tau) ...
    + gamma(1,2,0) / (lambda(1,0) - lambda(2,0)) * c_2i) .* exp(beta(1,tau));

corr2 =  @(tau) c_1i * gamma(2, 1, tau) ./ (lambda(1, tau) - lambda(2, tau)) .* exp(beta(1, tau)) - ...
    ( c_2i * quad(@(rho) gamma(1,2, rho) .* gamma(2,1,rho) ./ (lambda(1, rho) - lambda(2,rho)), 0, tau) ...
    + gamma(2,1,0) / (lambda(1,0) - lambda(2,0)) * c_1i) .* exp(beta(2,tau));

tau = linspace(0, T, 30);
num_i = deval(num, tau / delta);

c = zeros(size(A(0),1), length(tau));
c_c = c;
c1_e = zeros(1, length(tau));
c2_e = c1_e;
c1_corr_e = c1_e;
c2_corr_e = c2_e;
c1_num = c1_e;
c2_num = c1_e;
for i = 1:length(tau)
    
    c1_e(i) = c1(tau(i));
    c2_e(i) = c2(tau(i));
    c1_corr_e(i) = corr1(tau(i));
    c2_corr_e(i) = corr2(tau(i));
    c1_num(i) = adj_evec(1, tau(i))' * num_i(:,i);
    c2_num(i) = adj_evec(2, tau(i))' * num_i(:,i);
    
    
    c(:,i) = c1_e(i) * evec(1, tau(i)) + c2_e(i) * evec(2,tau(i));
    c_c(:,i) = (c1_e(i) + delta * c1_corr_e(i)) * evec(1, tau(i)) ...
        + (c2_e(i) + delta * c2_corr_e(i)) * evec(2,tau(i));
end

% plot the actual solution
subplot(2,1,1);
plot(num.x * delta, num.y(1,:), 'k', tau, c(1,:), 'b*');
hold on;
plot(tau, c_c(1,:), 'ro');
xlabel('$\tau$','interpreter','latex', 'fontsize',12);
ylabel('$\hat{c}_1$','interpreter','latex', 'fontsize',12);
subplot(2,1,2);
plot(num.x * delta, num.y(2,:), 'k', tau, c(2,:), 'b*');
hold on;
plot(tau, c_c(2,:), 'ro');
xlabel('$\tau$','interpreter','latex', 'fontsize',12);
ylabel('$\hat{c}_2$','interpreter','latex', 'fontsize',12);
l = legend('numeric','leading order','leading plus correction');
set(l,'interpreter','latex','location','best');


% plot the 'Fourier' coefficients 
figure;
subplot(2,1,1);
plot(tau, c1_num, 'k', tau, c1_e, 'b*', tau, c1_e + delta * c1_corr_e, 'ro');
subplot(2,1,2);
plot(tau, c2_num, 'k', tau, c2_e, 'b*', tau, c2_e + delta * c2_corr_e, 'ro');



% plot the relative error of the first component
num_i = deval(num, tau / delta, 1);
figure;
semilogy(tau, abs((num_i - c(1,:)) ./ num_i), 'b');
hold on;
semilogy(tau, abs((num_i-c_c(1,:)) ./ num_i), 'r');
xlabel('$\tau$','interpreter','latex', 'fontsize',12);
ylabel('Relative error','interpreter','latex', 'fontsize',12);
l = legend('leading order','leading plus correction');
set(l, 'interpreter','latex','location','best');


%--------------------------------------------------------------------------
% corrects the indexing... (index = 1 is largest eigenvalue...)
function v = ind(i)
global A
v = size(A(0),1) - i + 1;

%--------------------------------------------------------------------------
% differentiation in time operator
function v = dt(object, tau)

v = (object(tau + 1e-5) - object(tau - 1e-5)) / 2e-5;

%--------------------------------------------------------------------------
% computes eigenvalue
function v = lambda(j, tau)
global A

for i = 1:length(tau)
    ev = eigs(A(tau(i)));
    v(i) = ev(ind(j));
end

%--------------------------------------------------------------------------
% computes eigenvector
function v = evec(i, tau)
global A

[ef, ev] = eigs(A(tau));
v = ef(:,ind(i));

%--------------------------------------------------------------------------
% computes adj eigenvector
function v = adj_evec(i, tau)
global A

[ef, ev] = eigs(A(tau)');
v = ef(:,ind(i)) / (ef(:,ind(i))' * evec(i,tau));


%--------------------------------------------------------------------------
% computes gamma_{j,k}
function v = gamma(j, k, tau)

for i = 1:length(tau)
    v(i) = -(adj_evec(j,tau(i)))' * dt(@(t) evec(k, t), tau(i));
end

%--------------------------------------------------------------------------
% computes beta_j
function v = beta(j, tau)
global delta

for i = 1:length(tau)
    v(i) = quad(@(T) lambda(j, T) / delta + gamma(j, j, T), 0, tau(i), 1e-10);
end

